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AMSRD-ARL-RO-OI  Proposal  Number:  45508-EG 


J.  N.  Reddy 

Department  of  Mechanical  Engineering 
Texas  A&M  University 


EXECUTIVE  SUMMARY 

Most  structural  components  encountered  in  anny  vehicles  and  annor  can  be  classified  as  beams, 
plates,  or  shells  for  analysis  purposes.  While  these  structural  elements  are  designed  to  function 
properly  under  thermo-mechanical  loads  encountered  in  their  use,  they  do  develop  high  stresses 
and  experience  high  vibration  frequencies  that  may  make  them  non-functional  in  actual  service 
conditions.  The  objective  of  this  research  is  to  develop  consistent  plate  and  shell  theories  and 
associated  computational  framework  for  linear  and  non-linear  problems  of  structural  dynamics  in 
which  localized  high  gradients  of  the  solutions  are  resolved  accurately  and  time  accuracy  of  the 
solution  is  assured  at  all  stages  during  the  evolution.  Crucial  importance  of  this  framework  will 
be  demonstrated  computationally  through  well  known  benchmark  model  problems  in  the  area  of 
solid  mechanics  with  special  focus  on  composite  structures.  The  developed  methodology  and  the 
resulting  infrastructure  with  its  applications  to  solid  and  structural  mechanics  problems  should 
provide  highly  reliable,  robust  and  accurate  computational  technology  to  the  United  States  Army 
Laboratories.  The  specific  objectives  of  this  research  were: 

•  Develop  accurate  and  consistent  structural  theories  and  associated  finite  element  models 
of  plates  and  shells  that  account  for  transverse  shear  deformation  and  illustrate  the 
accuracy  using  benchmark  plate  and  shell  problems. 

•  Develop  mixed  and  least-squares  finite  element  models  of  the  refined  theories  for  the 
analysis  of  plates  and  shells. 

In  the  following  pages  a  technical  discussion  of  the  scientific  progress  made  and 
accomplishments  are  summarized  in  two  parts: 

1 .  A  robust  shell  finite  element  for  nonlinear  analysis  of  composite  and  functionally  graded 
shells,  and 

2.  Mixed  least-squares  finite  element  models  for  bending  and  vibration  of  plates. 


2 


1.  LARGE  DEFORMATION  ANALYSIS  OF  SHELLS 


1.1  Introduction 

Composite  shells  have  been  of  great  interest  in  many  engineering  applications.  Composites 
made  up  of  fiber-reinforced  laminae  that  are  bonded  together  (Reddy,  2004)  are  particularly 
attractive.  A  typical  lamina  is  often  characterized  as  orthotropic  with  the  principal  material 
directions  of  each  lamina  coinciding  with  the  fiber  direction  and  transverse  to  it.  As  required  in  a 
design,  by  changing  the  material  type,  fiber  orientation,  or  thickness,  the  designer  can  tailor  the 
different  properties  of  a  laminate  to  suit  a  particular  application.  Despite  their  multiples 
advantages,  laminated  composites  exhibit  a  serious  shortcoming  due  to  concentrations  of 
stresses,  as  well  as  in-surface  displacements,  caused  by  the  piece-wise  variation  of  the  material 
properties  through  the  thickness  of  the  shell.  Consequently,  a  special  material  named 
“functionally  graded  materials”  (FGMs)  was  proposed  by  Koizumi  (1997)  and  Yamanouchi  et  al. 
(1990),  in  which  the  material  properties  vary  smoothly  and  continuously  from  one  surface  to  the 
other.  These  materials  are  inhomogeneous  and  made  from  isotropic  components.  The  gradation 
of  the  material  properties  through  the  thickness  avoids  jumps  or  abrupt  changes  on  the  stress  and 
displacement  distributions  of  any  thin-walled  structure. 

In  some  applications  shell  structures  can  experience  large  elastic  deformations  and  finite 
rotations.  Geometric  nonlinearity  plays  an  essential  role  in  the  behavior  of  the  shell,  especially 
when  it  reaches  large  deformations.  Previous  studies  show  that  laminated  shells  exhibit  drastic 
changes  in  their  bending  response  (Ba§ar  et  al.,  1993;  Vu-Quoc  and  Tan,  2003;  Balah  and  Al- 
Ghamedy,  2002).  Even  for  homogeneous  and  isotropic  shells  we  observe  an  unpredictable 
behavior  (Simo  et  al.,  1990;  Sansour  and  Kollmann,  2000).  Therefore,  it  is  of  vital  importance  to 
study  the  nonlinear  response  of  potentially  inhomogeneous  materials  such  as  functionally  graded 
shells. 

This  paper  is  motivated  by  the  lack  of  studies  found  in  the  literature  that  addresses  large 
deformation  analysis  for  FGM  shells.  A  review  of  technical  articles  shows  that  few  studies  have 
been  carried  out  to  investigate  the  nonlinear  bending  response  of  plates  and  shells.  Most  of  them 
use  von  Kannan  or  Sanders  theories  which  are  restricted  to  moderately  small  deformations.  We 
cite  the  papers  of  Na  and  Kim  (2005),  who  examined  the  effect  of  thennal  loading  and  uniform 
pressure  on  the  bending  response  of  FGM  plates;  and  Yang  and  Shen  (2003a, b),  who  analyzed 
the  nonlinear  bending  and  postbuckling  behavior  for  FGM  plates  under  thermomechanical  load 
with  various  boundary  conditions.  Woo  and  Meghid  (2001)  provided  an  analytical  solution  for 
large  deflection  FGM  plates  and  shells  under  mechanical  and  thermal  loading;  while  Ma  and 
Wang  (2003)  examined  the  axisymmetric  large  deflection  bending  and  thermal  postbuckling  of 
FGM  circular  plates  subjected  to  mechanical  and  thermal  loading.  Both  articles  are  based  on  the 
classical  von  Karman  plate  theory. 

Moreover,  Reddy  and  Chin  (1998)  analyzed  the  dynamic  thennoelastic  response  of 
functionally  graded  cylinders  and  plates.  Praveen  and  Reddy  (1998)  carried  out  a  nonlinear 
thennoelastic  analysis  of  functionally  graded  ceramic-metal  plates  using  a  finite  element  model 
based  on  the  FSDT.  Thermomechanical  buckling,  as  well  as  bending  and  free  vibration  analysis, 
of  FGM  plates  can  be  found  in  the  articles  by  Reddy  and  Arciniega  (2006a, b).  Further  studies  of 
bending  and  vibration  analyses  of  FGMs  plates  can  be  found  in  the  articles  of  Reddy  (2000),  and 
Della  Croce  and  Venini  (2004). 
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On  the  subject  of  computational  models  for  shell  structures,  we  focus  our  attention  on  tensor- 
based  finite  element  models  (Harte  and  Eckstein,  1986).  This  approach  is  able  to  detennine  all 
properties  of  the  shell’s  differential  geometry  exactly.  Additional  errors,  introduced  by 
approximating  the  geometry  of  the  midsurface  of  the  shell  (as  in  continuum-based  finite  element 
models),  are  prevented  from  the  beginning.  Previous  works  of  the  authors  using  tensor-oriented 
finite  element  formulations  for  linear  analysis  of  laminated  shells  can  be  found  in  Arciniega  and 
Reddy  (2005),  and  Reddy  and  Arciniega  (2004). 

In  this  paper,  a  large  deformation  analysis  for  functionally  graded  shells  is  presented.  The 
formulation  is  based  on  the  first-order  shear  deformation  theory  with  seven  independent 
parameters  (Sansour,  1995;  Bischoff  and  Ramm,  1997)  where  no  plane  stress  assumption  is 
required  (3D  constitutive  equations).  A  tensor-based  finite  element  model  is  developed  using 
high-order  Lagrange  elements  to  preclude  membrane  shear  locking.  The  gradation  of  the 
material  properties  of  the  FGM  shell  is  considered  through  the  thickness.  The  material  stiffness 
tensor  is  obtained  by  Gauss  integration.  Numerical  results  are  presented  for  typical  benchmark 
problems  with  applications  to  functionally  graded  shells. 

1.2.  Theoretical  Formulation 

The  shell  theory  will  be  briefly  discussed  here.  For  a  detailed  development,  one  can  consult  the 
paper  of  Arciniega  and  Reddy  (2006)  and  references  herein.  The  mathematical  background 
utilized  in  the  following  derivation  is  given  in  the  books  of  Naghdi  (1963)  and  (1972),  Green  and 
Zerna  (1968),  and  Pietraszkiewicz  (1979). 

Let  us  introduce  in  the  region  B  fi(B  r)  a  converted  curvilinear  coordinate  system 
{ O' } ,  i  =  1,2,3,  such  that  the  surface  O  '  =  0  defines  the  midsurface  M  A,(M  ()  of  the  region 
B«(B  ,).  The  coordinate  O'  is  the  measure  of  the  distance  between  points  P  G  B*  (Pe B,) 
and  M  G  M  R  ( M  gM  (),  with  -  h!2<6 3  <h/2,  where  h  is  the  thickness  of  the  shell  (Fig.  1). 


Fig.  1.  Shell  continuum  in  the  reference  and  current  configurations. 
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Consider  the  motion  %(X,t)  of  the  shell  body  B  from  the  reference  configuration  B  R  to 
the  current  configuration  B, .  Since  a  convected  coordinate  system  {& }  has  been  adopted, 
geometric  quantities  of  the  region  B  t  are  analogous  to  those  defined  in  B  R .  In  the  Lagrangian 
description,  the  displacement  of  the  particle  X  from  the  reference  configuration  to  the  current 
configuration  is  given  by  the  vector  v(X,t) ,  i.e. 


(1) 


v(X,t)  =  Z(X,t)-X  =  x-X 

=  v‘g  l  =  Vjgj 

wherein  the  last  line  is  in  component  form  with  respect  to  the  region  B  R . 

We  introduce  the  first  kinematical  assumption  for  the  shell  model:  “the  displacement  field  is 
considered  as  a  linear  expansion  of  the  thickness  coordinate  around  the  midsurface.  The 
transverse  displacement  is  parabolic  through  the  thickness  of  the  shell”. 

This  assumption  implies  that 


where 


v(  ea )  =  u(  ea ) + e\{ea )  +(#3)2  y(  ea ) 


u(0a)  =  uiai,  (p(0a)  =  <piai,  \f{0a)  =  y/^ 


(2) 


(3) 


The  underlined  tenn  of  equation  (2)  is  included  to  avoid  Poisson  locking  (Bischoff  and  Ramm, 
1997). 

The  position  vector  of  the  deformed  shell  can  be  obtained  substituting  equation  (2)  into  (1). 
Thus 

x=r+0%+(03)2\\r  (4) 

where  r=r+u  and  a3  =a3+  cp.  The  vector  tp  is  also  called  difference  vector  (change  of  the 
director  of  the  midsurface).  The  director  a3  is,  in  general,  neither  a  unit  vector  nor  orthogonal  to 
M  t .  The  configuration  of  the  shell  is  uniquely  detennined  by  the  displacement  vector  u  of  the 
midsurface  together  with  the  difference  vector  (p  and  the  additional  variable  tj/,  or  by  seven 

independent  components  of  these  vectors  (Sansour,  1995). 

We  now  introduce  the  Green  strain  tensor  E  as  a  measure  of  the  strain  for  a  material 
description 


E  =  I(C  — G) 


(5) 


where  C  =  FrF  is  the  right  Cauchy-Green  tensor,  G  =  gif  g'  ®  g7  is  the  Riemannian  metric  in 
the  reference  configuration  and  F  =  g,  0  g'  is  the  defonnation  gradient.  We  define  the  covariant 
space  and  surface  base  vectors  in  the  current  configuration  as  g,  and  a) ,  respectively. 

The  shifter  tensor  p  is  a  two-point  tensor  which  relates  the  region  B  A,  to  the  reference 

midsurface  M  R  and  it  is  useful  to  define  the  tensor  E  as 
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(6) 


E  =  0*  (E)  =  nrEn 

where  <f>*(°)  is  the  pull-back  operator. 

The  tensor  E  can  be  expanded  as  a  function  of  the  thickness  coordinate,  i.e. 

E  =  £°  +g3E1+(fl3)2E2  +  (fl3)V  +  (fl3)V  (7) 

The  second  assumption  for  the  shell  model  asserts  that:  “quadratic  and  higher-order  terms  of 
E  ,  underlined  in  equation  (7),  are  negligible”.  Then,  we  arrive  to  the  following  decomposition 

e°  =  a“  <g>  (a"  ®  a3  +a3  ®  a")  +  ®  a3 

e1  =  £(lp  a“  ®  a^  +  £™  (a“  ®  a3  +a3  <g>  a“  )  +  eQ  a3  ®  a3 

where  £„},£„ \  and  £p  are  functions  of  the  triple  fu,  (p,  \)/j .  After  some  manipulations  we  can 
write  them  in  terms  of  the  seven  components  of  the  displacement  field  (Habip,  1965),  i.e. 

£aj  =  \{Ua\p  +  UP\ a  ~  2KpU 3  +  ~b}U 3UA\a  ~KU3UA\p 

+CafMzf  +UXaUXP  +baUlUXP  +  bpUAUXa  +  babpUAUy) 

£a}  =  ~^Pa\p  ~>r(Pp\ a  ~2bap(p3  ~~bpUA\a  ~baUA\p  +2cctPU 3  +  a  1'UA\a<Py\p 
WruAp(pAa  -bp(p,uAa  -bfau^  -bpU,(pMa  —bxau%(px^  +  2 

+U3,a(Pi,p  +uxp(Pxa  +  ba<PAuxp  +bp(pxuXa  +  b^(Ux  (pifj  +  b'pUx(pia 

+bablUA(Py+bpKUA(Py) 

£a3  =  2  (^“  +  W3,a  +  baUA  +  aA]'UA\a(Py  ~  ba<PAU 3  +  ^3  W3,a  +  ^3  ) 

<3  =  |(^3.«  +  +  ^3  +  2^3M3,«  +  2^X  ) 

40)  =  ^(2^3  +  aXr(p,(pr  +  (p3)2 ) 

4  =2(^3  +^3^3)  (9) 

where  cap  =  baJjp  is  the  covariant  third  fundamental  form  of  the  reference  surface.  Note  that  the 
component  £p  vanishes  when  !//,  =  0  (6-parameter  formulation). 

The  second  Piola-Kirchhoff  stress  tensor  is  used  for  the  Lagrangian  formulation  and  is 
energetically-conjugate  to  the  rate  of  Green  strain  tensor  E  (Reddy,  2004).  Like  E  ,  the  second 
Piola-Kirchhoff  stress  tensor  S  is  transformed  to  the  midsurface  M  R  by 

S  =  p-1Sp”r  =0*(S)  (10) 

which  is  the  pull-back  operator  of  the  contravariant  tensor  S  . 
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Let  M"  denote  the  stress  resultant  tensor  which  is  a  symmetric  tensor.  The  tensor  M"  is 
defined  as 


[m0,M']  =  J/!/2Ji,6»3]s  ludd*  (11) 

The  scalar  quantity  //  is  the  detenninant  of  the  shifter  tensor  p .  The  stress  resultant  tensors  are 

also  energetically-conjugate  to  the  strain  resultants  s' .  The  stress  resultant  tensors  may  be 
decomposed  in  component  fonn  as 

(0)  (0)  (0) 

M°  =  NaP  aa  <g)  +  Q "3  (aa  <g  a3  +a3  <g  aa)  +T33  a3  <g  a3 

,  (1)  «  (D  ,  0),-  n  2) 

M  =  Nap  aa®ap  +  Q  (aa  <g  a3  +a3  <g  aa)  +T: a3  <g>  a3 

(«)  „  («)  («) 

where  Nap,Qa~  and  T  are  membrane,  shear  and  stretching  components,  respectively. 

1.3.  Functionally  Graded  Shells 

In  this  section  we  consider  a  hyperelastic  and  inhomogeneous  shell.  The  shell  structure  can 
undergo  large  defonnations  (rotations  and  displacements)  while  the  material  response  remains  in 
the  elastic  regime.  We  also  consider  the  relation  between  the  second  Piola-Kirchhoff  stress 
tensor  S  and  the  Green  strain  tensor  E  is  linear.  It  implies  that 

S  =  C-E  (13) 

where  C  is  the  fourth-order  elasticity  tensor.  The  tensor  C  is  represented  in  convected 
coordinates  as 


C  =  CiJk!  g,.  ®  gj  ®  gk  ®  g,  (14) 

where  the  components  of  C  satisfy  the  following  symmetry  conditions 

Cijkl  =  Cjikl  =  Cij,k  =  CkHj  ( 1 5) 

Functionally  graded  materials  (FGMs)  are  a  special  kind  of  composites  in  which  the  material 
properties  vary  smoothly  and  continuously  from  one  surface  to  the  other.  These  materials  are 
microscopically  inhomogeneous  and  are  typically  made  from  isotropic  components.  One  of  the 
main  advantages  of  FGMs  is  that  it  mitigates  severe  stress  concentrations  and  singularities  at 
intersections  between  interfaces  usually  presented  in  laminate  composites  due  to  their  abrupt 
transitions  in  material  compositions  and  properties.  Applications  of  FGMs  are  extensive 
especially  in  high-temperature  environments  such  as  nuclear  reactors,  chemical  plants  and  high¬ 
speed  spacecrafts. 

The  materials  in  the  bottom  and  top  surfaces  are  usually  metal  and  ceramic  respectively  (Fig. 
2).  Material  properties  at  a  point  X  are  given  by  a  combination  between  metal  and  ceramic 
constituents,  i.e.  by  the  weighted  average  of  the  moduli  of  the  constituents,  namely 

Cu{Oi)  =  CUcfc+CUmfm  (16) 

where  the  subscripts  m  and  c  refer  to  the  metal  and  ceramic  constituents  and /is  the  volume 
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fraction  of  the  phase.  The  symbol  U5  denotes  a  generic  material  property. 


Fig.  2.  Functionally  graded  shell. 


The  volume  fractions  of  the  ceramic  fc  and  metal  fm  corresponding  to  the  power  law  are 
expressed  as  (Reddy,  2000;  Praveen  and  Reddy,  1998;  Reddy  and  Chin,  1998) 


fc  = 


Z  1 

7i  +  2 


f  =1  -f 

J  m  J  c 


(17) 


where  n  is  the  volume  fraction  exponent  which  takes  values  greater  than  or  equal  to  zero.  The 
value  of  n  equal  to  zero  represents  a  fully  ceramic  shell.  Conversely,  we  have  a  fully  metal  shell 
as  n  tends  to  infinity  (Fig.  3). 


6>V  h 


Fig.  3.  Variation  of  the  volume  fraction  function  fc  through  the  dimensionless  thickness  for 
different  values  of  power-law  index  n . 
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The  components  of  the  elasticity  tensor  Cijkl  ( O' )  are  functions  of  the  thickness  coordinate. 
They  can  be  written  in  terms  of  the  convected  base  vectors  as 


C  =  Cijkl  03)g;®g/®g,®g/ 


(18) 


which  can  be  arranged  in  a  matrix  [C 


ijkl 


[[6x6 


such  that 


[cm\  = 


1111 


1122 


1133 


c 

c 

c 


1122 


2222 


2233 


c 

c 

c 


1133 


2233 


3333 


0 

0 

0 


0 

0 

0 


0 

0 

0 


c 


0 

0 

0 

2323 

0 

0 


c 


0 

0 

0 

0 

1313 


0  c 


0 

0 

0 

0 

0 

1212 


(19) 


6x6 


The  components  Cijkl  at  each  0  are  functions  of  only  two  independent  variables,  then 


Cim  =  C  LILA  _  _ 


~i2222 


3333 


^~rll22  _  ^~fll 33  ^~iAA5i  _ 


-fl  133 


2233 


£tl212  _  £1313  _  ^2323  _ 


i-v) 
(l  +  v)(l-2v) 
E(0^)v 


(l  +  v)(l  — 2v) 

e(o^ 


(20) 


2(l  +  v) 


where  E  ( 0  j  =  Ecfc  +  Emfm .  The  Poisson’s  ratio  v  is  considered  constant  through  the  thickness. 
Hence 


_  r'ijkl  r  sjijM 
^  cm  J  c  '  m 

where  Cf'  =  Cf* -Cf  and  fc,fm  are  given  in  (17). 


(21) 


1.4.  Weak  Formulation 


The  finite  element  framework  is  based  on  the  principle  of  virtual  work.  Our  analysis  is  restricted 
to  static  cases.  The  virtual  work  statement  is  nothing  but  the  weak  form  of  the  equilibrium 
equations  and  it  is  valid  for  linear  and  nonlinear  stress-strain  relations  (Reddy,  2002). 

The  abstract  configuration  solution  of  the  shell  is  denoted  by  the  set 


C  ={<D  =  (u,(p,\j/) 


0:A  e 


xiirxj 


(22) 


where  A  is  the  parametric  space  of  the  midsurface.  Note  that  OgC  contains  the  same 
amount  of  three-dimensional  infonnation  as  Eq.  (2)  to  locate  at  any  time  arbitrary  points  in  the 
three-dimensional  shell. 
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We  express  the  weak  formulation  as 


where 


G  (o,  <5o)  =  gmt  (o,  <5o) + gext  (o,  <5o)  =  o 

£int(0,<50)  =  f  (M°-  Sz0  +M1-  5zl)d£l 
Jm  r 

£ext(0,<50)  =  -JM  (p<5u  +  l<5cp  +  k<5\j/)</Q 

—  f  (p'-<5u  +  r'-<5cp  +  ks-<5\|/)<A 

J  cM  »  — 


(23) 

(24) 

(25) 


For  hyperelastic  materials,  the  static  part  of  the  weak  form  of  the  equilibrium  equations  is  the 
first  variation  of  an  elastic  potential  energy  function.  This  statement  is  known  as  the  principle  of 
minimum  total  potential  energy  (Reddy,  2002).  We  define  the  elastic  potential  function 

Id  ( °) :  C  — ►  R  as 


n(0)=  f  p^dV—  f  (p-u +  l-(p  +  k-\|/)</Q 

Jbs  Jm  „  — 

-f  (ps u  +  ls -(p  +  k*  y)ds 

6 M  R  — 

The  first  variation  of  the  potential  energy  is  given  by 


(26) 


g  (o,  <5o)  -  sn  ( o,  <5o)  -  mi  (o)  [<5o  ]  =  o  (27) 

To  solve  the  nonlinear  equations  is  to  use  the  incremental/iterative  method  of  Newton- 
Raphson.  This  procedure  requires  a  linearization  of  the  weak  form  generating  recurrence  update 
formulas.  The  linearization  process  relies  on  the  concept  of  directional  derivatives  (Hughes  and 
Pister,  1978;  Bonet  and  Wood,  1997).  We  assume  that  the  external  forces  are  conservative 
(independent  of  O  ).  Applying  that  procedure  to  equation  (23)  we  obtain 


eg  (O,  <50;  AO)  =  g  (O,  <50 )  +Dg(  O,  <50 )  [AO]  +  o  ( AO)  (28) 

where  the  underlined  term  is  called  consistent  tangent  operator.  Furthermore,  we  can  write  the 
tangent  operator  as 


£>£(0,<50)[A0]  =  V£(0,<50)-A0  (29) 

since  <50  remains  constant  during  the  increment  AO  . 

The  iterative  solution  procedure  goes  as  follows:  given  a  configuration  Oa  gC  , 
corresponding  to  iteration  k,  solve  the  linearized  system 

£(0a,<50)  +  W(0a,<50)-A0a  =0  (30) 

where  AO'  is  the  incremental  change  in  the  configuration  of  the  shell.  This  increment  is  used  to 
update  the  shell  configuration  O  — ►  Oa+1  gC  .  Namely 

Oa+1=Oa+AOa  (31) 
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Notice  that  the  use  of  the  triple  (u,cp,\|/)  preserves  the  additive  structure  of  the  configuration 
update  of  the  shell. 

The  consistent  tangent  operator  is  decomposed  in  two  parts:  the  material  tangent  operator  and 
the  geometric  tangent  operator.  Thus 

D  g  (O,  ^O)  [AO]  =  Dmg  (O,  AO)[AO]  +  DgQ  (O,  <SO)  [AO]  (32) 

The  contribution  of  the  external  forces  vanishes  because  they  are  conservative.  The  first  term 
which  is  the  material  part  is  given  by 

DmQ  (0,£0)  [AO]  =  f  V(DM"[AO]-<&")rfQ  (33) 

M  p  „ 

R  n=0 

and  the  geometric  part  by 

DaQ  (O,  (50)  [AO]  =  f  V(M"-jD^£"[AO])JQ  (34) 

The  material  part  of  the  tangent  operator  results  from  the  directional  derivative  of  the  stress 
resultants.  After  some  manipulations  we  obtain 

Z>M''(0)[A0]  =  f"''2jU(03y+JC-AeJ  dO3  (35) 

where  C  is  the  pull-back  of  the  contravariant  fourth-order  elasticity  tensor  C .  Substituting  (35) 
into  (33)  we  arrive  to 


Dm  Q  (O,  (50)  [AO]  =  f  VV(^‘-B(I^-A  sJ)dQ,  (36) 

/—O  y=0 

where  Aej  is  can  be  easily  calculated.  The  components  of  the  fourth-order  tensor  B(i)  are  the 
material  stiffness  coefficients  of  the  shell  and  are  defined  as 


B(t)  =  fh'2  jU(03)kCd03,  A:  =  0,1,2  (37) 

J  —h/2 

and  are  computed  by  Gauss  integration. 

The  computation  of  virtual  internal  energy  gint  and  the  tangent  operator  is  not  a  trivial  task. 

Even  for  isotropic  materials  these  expressions  have  an  extremely  complex  form  when 
displacements  and  rotations  are  large. 

Next,  the  finite  element  equations  are  obtained  by  interpolating  the  covariant  components  of 
the  kinematic  variables  in  terms  of  the  base  vectors  a" .  Namely 


u^(0)=  a\  (p/,p(0)=  'E<plU)NU)(£,rl) 


j= 1 


Wn\  _ 


¥  6 


)]](/a!,i.v'/’(c.//) 


j= 1 


(38) 
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where  (u^f  ,^30))  denote  the  nodal  values  of  the  kinematic  variables. 

We  then  arrive  to  a  system  of  highly  nonlinear  algebraic  equations  which  can  be  written  in 
matrix  form  by  means  of  the  stiffness  and  tangent  matrices.  The  solution  is  carried  out  by 
subroutines  written  in  FORTRAN. 

1.5.  Numerical  Examples 

In  this  section,  numerical  results  obtained  by  the  model  developed  herein  are  presented  for  shell 
structures.  Typical  benchmark  problems  for  isotropic  and  homogeneous  shells  are  investigated 
for  bending  behavior  of  their  counterparts  functionally  graded  shells. 

Regular  meshes  of  Q25,  Q4 9  and  QS 1  high-order  elements  with  seven  degrees  of  freedom 
per  node  were  utilized  in  the  finite  element  analysis  (see  Table  1).  By  increasing  the  p  level  or 
refining  the  finite  element  mesh,  we  mitigate  locking  problems.  Full  Gauss  integration  rule  is 
employed  in  all  examples. 


Table  1.  Number  of  degrees  of  freedom  per  element  for  different  p  levels 


Element 

p  level 

FSDT  (DOF) 

Q4 

1 

28 

Q9 

2 

63 

Q25 

4 

175 

Q49 

6 

343 

Q81 

8 

567 

Roll-up  of  a  functionally  graded  plate  strip 

We  consider  a  FGM  plate  strip  subjected  to  a  bending  distributed  moment  on  the  other  end  (Fig. 
4).  The  isotropic  and  homogeneous  counterpart  has  been  considered  Simo  et  al.  (1990)  as  well  as 
Betsch  et  al.  (1998).  This  problem  is  good  to  test  the  capability  of  the  finite  element  model  to 
simulate  large  rotations  on  shells. 


Fig.  4.  Cantilever  FGM  plate  strip  under  end  bending  moment. 
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The  material  properties  and  geometry  of  the  plate  are 


Em  =  0.7  xlO9,  Ec  =1.51xl09,  v  =  0.3 
L  =  12.0,  6  =  1.0,  A  =  0.1 
=65886.17926 

Figures  5  and  6  depict  tip  displacements  of  the  cantilever  strip  plate  versus  the  end  bending 
moment  for  various  volume  fraction  exponents  n  (from  fully  ceramic  to  fully  metal).  We  utilize  a 
regular  mesh  of  1x8(925  elements  for  the  finite  element  discretization.  The  Newton  method 
exhibits  a  good  rate  of  convergence  until  some  displacement  level  and  then  it  diverges  (for 
inhomogeneous  shell  cases).  It  is  not  clear  for  the  authors  why  this  problem  happens.  It  seems 
that  for  these  cases  we  do  not  have  real  solutions.  However,  before  arriving  to  any  conclusion 
further  studies  are  needed. 


Fig.  5.  Tip-deflection  —  uKl>  vs.  end  moment  M for  the  FGM  plate  strip. 


Figure  7  shows  the  undeformed  and  deformed  configuration  of  a  FGM  strip  plate  for  various 
load  stages  and  «  =  1.0.  The  plate  shows  large  rotations  beyond  180°  with  deformed 
configurations  similar  to  the  homogeneous  case. 

Annular  FGM  plate  under  end  shear  force 

We  analyze  an  annular  FGM  plate  subjected  to  a  distributed  transverse  shear  force  (Fig.  8).  This 
benchmark  problem  was  considered  for  homogeneous  and  isotropic  plates  by  Buchter  and  Ramm 
(1992)  and  Sansour  and  Kollmann  (2000);  and  for  multilayered  composites  by  Arciniega  and 
Reddy  (2006).  The  material  properties  are  the  same  as  the  last  example  and  will  be  used  in  all 
examples.  The  geometric  quantities  are  given  by 
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R.  =  6,  Re  =10 ,h  =  0.03 
for  a  maximum  distributed  force  of  qmax  =  20.0 . 


Fig.  6.  Tip-deflection  u<3>  vs.  end  moment  M  for  the  FGM  plate  strip. 


Fig.  7.  Deformed  configurations  of  a  FGM  plate  (n  =  1.0)  under  end  bending  moment  (load 
values  M/Mref  =  0.075, 0. 1 5, . . . ,  0.6, 0.625 ). 
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Fig.  8.  Annular  FGM  plate  strip  under  transverse  end  shear  force. 

The  plate  is  modeled  by  polar  coordinates.  A  regular  mesh  of  1x5049  elements  (p  level 
equal  to  6)  is  used  in  the  present  analysis.  Computation  is  performed  by  the  Newton-Raphson 
method  with  80  load  steps  and  convergence  tolerance  for  the  residual  forces  of  1.0x10  4 . 

The  shear  load  versus  displacement  curves  for  two  characteristic  points  are  depicted  in 
Figures  9  and  10.  The  deformed  configurations  of  a  FGM  annular  plate  for  various  load  levels 
and  n  =  2.0  is  shown  in  Fig.  1 1.  It  is  clear  that  the  plate  undergoes  large  displacements  at  the 
corresponding  loading  of  F  —  80  . 


Fig.  9.  Transverse  displacement  curves  at  point  A  vs.  shear  force  F  =  4q  of  the  cantilever 
annular  plate  strip. 
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Fig.  10.  Transverse  displacement  curves  at  point  B  vs.  shear  force  F  =  4q  of  the  cantilever 
annular  plate  strip. 


Fig.  11.  Deformed  configurations  of  a  FGM  plate  strip  (n  =  2.0)  under  transverse  end  shear 
force  (load  values  F  =  10, 20,..., 80). 
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Pull-out  of  a  functionally  graded  cylindrical  shell 

The  functionally  graded  cylindrical  shell  with  free  ends  is  subjected  to  two  opposite  loads  (Fig. 
12).  The  homogeneous  case  was  considered  by  Brank  et  al.  (1995)  and  Sansour  and  Kollmann 
(2000),  among  others.  The  following  geometrical  data  is  used  in  the  analysis 

L  =  10.35,  R  =  4.953,  h  =  0.094 

An  octant  of  the  shell  is  modeled  using  2x2  (981  elements  which  is  enough  to  overcome 
locking  problems.  The  Newton-Raphson  method  with  80  load  steps  is  utilized  with  equal  load 
steps  of  60000.  The  adopted  error  tolerance  for  the  residual  was  1  .Ox  1 0'5. 


c 


Fig.  12.  Pull-out  of  a  FGM  cylinder  with  free  edges. 

Figures  13  to  15  show  the  radial  displacements  at  points  A,  B  and  C  of  the  shell,  respectively. 
Convergence  rates  for  this  example  are  quite  good  (3  to  5  iterations  per  load  step).  As  expected, 
bending  response  of  FGM  cylinders  lies  in  between  of  the  fully  ceramic  and  fully  metal  shells. 
The  deformed  configurations  for  a  FGM  cylindrical  shell  is  depicted  in  Fig.  16  for  P  =  5.1xl06 
and  n  — 1.0. 


FGM  hemisphere  under  internal  pressure 

The  last  example  considered  is  a  cylindrical  FGM  shell  under  internal  pressure  (Fig.  17).  This  is 
not  a  following  loading  (independent  of  the  displacements). The  cylinder  has  fixed  boundary 
conditions  on  both  ends.  The  geometric  data  is  as  follows: 

a  =  20.0,  R  =  5.0,  h  =  0.01 
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0.0  0.5  1.0  1.5  2.0  2.5  3.0 

Radial  deflection  at  point  A 


Fig.  13.  Radial  displacements  at  point  A  (n<3>)  vs.  pulling  force  of  a  FGM  cylinder  with  free 
edges. 


Fig.  14.  Radial  displacements  at  point  B  (—it  ;3> )  vs.  pulling  force  of  a  FGM  cylinder  with  free 
edges. 
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0.0  0.5  1.0  1.5  2.0  2.5  3.0  3.5  4.0 

Radial  deflection  at  point  C 


Fig.  15.  Radial  displacements  at  point  C  (— «<3>)  vs.  pulling  force  of  a  FGM  cylinder  with  free 
edges. 


P  i 


Fig.  16.  Deformed  configurations  of  the  FGM  cylinder  under  pulling  forces.  Load  P  =  5.1xl06 
0  =  1.0). 
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FGM  hemisphere  under  internal  pressure 


The  last  example  considered  is  a  cylindrical  FGM  shell  under  internal  pressure  (Fig.  17).  The 
cylinder  has  fixed  boundary  conditions  on  both  ends.  The  geometric  data  is  as  follows: 
a  =  20.0,  R  =  5.0,  h  =  0.01 


C 


Fig.  17.  FGM  cylindrical  shell  under  internal  pressure. 


A  regular  mesh  of  2x208 1  elements  is  used  in  the  analysis.  We  take  advantage  of  the 
symmetry  of  the  shell  and  only  an  octant  of  the  shell  is  considered  as  the  computational  domain. 
Figure  18  shows  the  radial  deflections  at  the  central  point  versus  the  internal  pressure  for  FGM 
cylinders.  We  notice  that  FGM  cylinders  with  low  values  of  n  exhibit  stiffer  response  than  those 
with  high  volume  fraction  exponent  (more  metal  than  ceramic).  The  final  configuration  of  a 
FGM  cylinder  for  n  —  5.0  is  depicted  in  Fig.  19. 


0.0  0.5  1.0  1.5  2.0  2.5  3.0  3.5  4.0 

Radial  deflection  at  point  A 


Fig.  18.  Radial  deflection  at^4  vs.  pressure  load  (0  =  106g)  of  a  FGM  cylindrical  shell. 
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Fig.  19.  Deformed  configuration  of  a  FGM  cylindrical  shell.  Loading  q  =  3.6  xlO6  (n  =  5.0). 


1.6.  Conclusions 

In  this  paper  we  present  a  large  deformation  analysis  for  functionally  graded  shells.  We  consider 
a  through- the-thickness  variation  of  the  material  properties  of  the  FGM  shell  which  is  made  of 
two  isotropic  constituents.  The  gradation  of  properties  through  the  thickness  is  assumed  to  be  of 
the  power  law  type.  A  tensor-based  finite  element  model  is  developed  for  geometric  nonlinear 
analysis  of  the  shell.  This  approach  showed  to  be  reliable  and  efficient.  The  derived  first-order 
shell  theory  with  seven  parameters  with  exact  nonlinear  deformations  is  consistent  and  simple.  It 
incorporates  thickness  changes  in  the  model,  and  then  full  3D  constitutive  equations  are  utilized. 
A  family  of  High-order  Lagrangian  elements  was  introduced  to  avoid  membrane  and  shear 
locking  for  shells.  We  found  that  the  nonlinear  bending  response  of  FGM  shells  under 
mechanical  loading  lies  between  that  of  ceramic  and  metal  shells,  as  expected.  The  patterns  of 
load-displacement  equilibrium  curves  of  FGM  shells  are  similar  to  those  of  isotropic  and 
homogeneous  counterparts.  Numerical  examples  for  plates  and  cylindrical  shells,  presented 
herein,  illustrate  the  validity  of  the  present  approach  and  the  developed  formulation  for  FGM 
shells. 
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2.  MIXED  LEAST-SQUARES  FINITE  ELEMENT  MODELS  FOR 
STATIC  AND  FREE  VIBRATION  ANALYSIS  OF  PLATES 


2.1.  Introduction 

Finite  element  models  for  the  analysis  of  multilayered  composite  plate  and  shell  structures  have 
been  widely  developed  in  the  last  few  decades.  In  overview,  the  main  approaches  to  establish 
plate  and  shell  theories  have  been  in  the  framework  of  the  so  called  axiomatic  theories.  The 
formulations  differ  in  equivalent  single-layer  or  layerwise  variable  descriptions  and  also  in  the 
choice  for  the  unknown  variables,  resulting  in  displacement,  stress  or  mixed  fonnulations  [1-3]. 
Traditionally,  variational  principles  have  been  established  to  derive  governing  equations 
consistent  with  the  chosen  fonnulations.  The  widespread  displacement  fonnulations  usually 
relate  to  the  well-known  principle  of  virtual  displacement  and  the  alternative  mixed  formulations 
typically  derive  from  the  Flu-Washizu  or  the  Hellinger-Reissner  variational  principles  [4,5], 
Naturally,  classical  finite  element  models  that  were  originally  developed  for  one-layered 
isotropic  structures  were  extended  in  a  straightforward  manner  to  multilayered  plates  and  shells 
[6,7].  The  classical  lamination  theory  (CLT),  first-order  shear  defonnation  theory  (FSDT)  and 
high-order  theories  have  been  known  to  provide  a  sufficiently  accurate  description  of  the  global 
response  of  multilayered  structures,  as  long  as  thin  to  moderately  thick.  Understandably,  for  a 
detailed  response  of  individual  layers  or  local  phenomena  description  one  must  use  layerwise 
theories  (LWT)  or  the  so  called  zig-zag  theories  that  were  indeed  entirely  originated  and  devoted 
to  layered  structures.  In  fact,  one  crucial  issue  for  these  theories  is  the  fulfilment  of  the  Cz°- 
requirements.  Basically,  this  means  that  displacements  and  transverse  stresses  must  be  C°- 
continuous  functions  in  the  thickness  direction  due  to  interlaminar  compatibility  and  equilibrium 
reasons. 

The  motivation  for  the  proposed  finite  element  models  comes  from  earlier  works  on 
mixed  finite  element  fonnulations  based  on  least-squares  variational  principle.  Namely,  works 
by  Pontaza  and  Reddy  [8],  Pontaza  [9]  and  also,  Duan  and  Lin  [10].  Overall,  the  least-squares 
finite  element  formulations  have  shown  promising  theoretical  and  computational  advantages, 
both  in  fluid  and  in  solid  mechanics.  Specifically,  Pontaza  and  Reddy  [8]  developed  a  mixed 
model  based  on  least-squares  formulation  for  the  bending  of  single-layered  isotropic  plates, 
using  the  classical  plate  theory  and  first-order  shear  deformation  theory.  The  prospect  of 
extending  this  model  gave  rise  to  the  proposed  mixed  least-squares  FSDT  model  for  static 
analysis  of  laminated  composite  plates.  Then,  a  pioneer  attempt  to  use  least-squares  formulation 
in  modal  analysis  led  to  the  development  of  the  mixed  least-squares  FSDT  finite  element  model 
for  free  vibration  analysis  of  laminated  composite  plates. 

The  least-squares  fonnulations  as  any  weighted  residual  formulation  provide  an 
alternative  approach  to  the  weak  form  finite  element  models,  both  displacement-based  and 
mixed.  In  the  framework  of  FSDT  weak  form  models,  displacement  formulations  are  known  to 
encounter  computational  difficulties  when  modelling  thin  plates.  The  finite  elements  become 
excessively  stiff,  which  results  in  an  erroneous  underprediction  of  displacements  in  static 
analysis  or  else  in  a  severe  overprediction  of  frequencies  in  free  vibration  analysis.  This 
phenomenon  is  known  as  shear-locking.  In  essence,  it  is  due  to  the  inability  of  shear  deformable 
elements  to  accurately  model  bending  within  an  element  under  a  state  of  zero  transverse  shearing 
strain.  Higher-order  elements  experience  relatively  less  locking,  but  sometimes  at  the  expense  of 
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a  slower  convergence.  Usually  shear-locking  problems  are  only  avoided  by  numerical  integration 
techniques.  Another  possibility  is  mixed  formulations,  where  in  addition  to  generalized 
displacements  the  stress  resultants  are  also  used  as  independent  variables  [11,12],  Mixed  finite 
element  models  based  on  weak  formulations  need  however,  that  the  finite  element  approximation 
spaces  satisfy  a  so  called  Inf-Sup  condition,  in  order  to  be  consistent  models  [13].  This  fulfilment 
is  in  general  known  to  be  rather  difficult  to  prove  analytically.  Furthermore,  mixed  weak  form 
models  yield  symmetric  but  not  positive-definite  stiffness  matrices,  adding  numerical  complexity 
to  the  models.  Alternatively,  within  weighted  residual  fonnulations,  least-squares  finite  element 
models  are  distinctive  for  being  solely  based  on  the  idea  of  minimizing  the  error  introduced  in 
the  approximation  of  the  governing  equations.  Then,  the  benefit  of  using  least-squares  variational 
principle  along  with  mixed  formulations  is  that  it  leads  to  a  variational  unconstrained 
minimization  problem,  where  the  finite  element  approximation  spaces  can  be  chosen 
independently.  Therefore,  stability  requirements  such  as  the  Inf-Sup  condition  never  arise.  This 
is  precisely  the  theoretical  merit  of  mixed  least-squares  formulations  as  it  was  demonstrated  in 
the  aforementioned  works  on  this  matter  [8-10]. 

The  proposed  mixed  least-squares  finite  element  models  consider  the  FSDT  with 
generalized  displacements  and  stress  resultants  as  independent  variables,  using  equal-order 
interpolation,  for  either  static  or  free  vibration  analysis  of  laminated  composite  plates. 
Specifically,  high-order  C°  basis  functions  and  full  integration  are  used  to  develop  the  discrete 
finite  element  models,  since  it  was  established  to  be  the  appropriate  way  to  truly  minimize  the 
least-squares  functional.  In  fact,  Pontaza  and  Reddy  [8]  and  later  Pontaza  [9]  demonstrated  the 
exponential  decay  of  the  least-squares  functional  with  increasing  order  of  the  element. 
Furthermore,  the  mixed  least-squares  model  for  static  analysis  uses  the  classical  C°  Lagrange 
basis  functions,  whereas  the  model  for  free  vibration  analysis  developed  later  uses  instead  C° 
interpolant  polynomials  of  Gauss-Lobatto-Legendre  quadrature  points,  which  are  more  suitable 
basis  functions  for  high-order  elements  [14].  Both  mixed  least-squares  discrete  models,  once  the 
boundary  conditions  are  properly  imposed,  yield  a  symmetric  and  positive-definite  stiffness 
matrix.  This  is  another  benefit  of  mixed  least-squares  models  as  opposed  to  mixed  weak  form 
models,  which  is  computationally  preferable.  Most  interestingly,  the  pioneer  mixed  least-squares 
model  for  free  vibration  analysis  yields  a  quadratic  eigenvalue  problem  with  symmetric  matrices, 
which  is  rather  atypical  within  conservative  systems.  Ultimately,  both  proposed  models  exhibit 
excellent  predictive  capabilities  in  the  framework  of  the  FSDT  as  demonstrated  by  numerical 
examples  presented  hereafter.  In  particular,  it  is  also  shown  that  both  least-squares  models  based 
on  high-order  basis  functions  are  insensitive  to  shear-locking. 

This  report  is  outlined  as  follows.  It  starts  by  introducing  the  governing  equations 
consistent  with  the  mixed  FSDT  finite  element  models  for  both  static  and  free  vibration  analysis 
of  laminated  composite  plates.  Then,  the  proposed  models  are  derived  from  the  least-squares 
formulation  and  related  finite  element  specifics  are  addressed.  Selected  numerical  examples  are 
presented  to  assess  the  predictive  capabilities  of  the  mixed  least-squares  models  through  static 
and  free  vibration  analysis  of  four  laminated  composite  plates  with  different  boundary  conditions 
and  various  side-to-thickness  ratios.  Lastly,  the  overall  conclusions  are  discussed. 

2.2  Governing  Equations 

Consider  a  laminated  composite  plate  of  total  thickness  h  and  composed  of  N  orthotropic  layers, 
as  shown  in  Fig.  1.  Typically,  the  layers  are  unidirectional  fibre-reinforced  laminas  whose  in- 
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plane  material  coordinate  axes  are  parallel  and  transverse  to  the  fibres  direction.  Thus,  the 
orientation  of  the  Ml  layer  is  defined  by  an  angle  Ok  between  the  plate  coordinate  x  and  the  fibres 
direction.  In  the  xv-plane,  Q  represents  the  undeformed  midplane  of  the  plate  and  Q  the  open 
bounded  region  with  the  boundary  T-dQ.  The  z-axis  is  taken  positive  upward  from  the 
midplane.  Specifically,  the  Mi  layer  is  located  between  the  interfaces  z  =  zk  and  z  =  zk+l  in  the 
thickness  direction. 


Figure  1.  Notations  for  a  laminated  composite  plate. 

As  previously  mentioned,  the  adopted  mixed  formulation  uses  the  generalized  displacements  and 
stress  resultants  as  independent  variables.  Accordingly,  for  static  analysis  of  laminated 
composite  plates  under  a  transverse  load  q(x,y),  the  governing  equations  by  the  FSDT  are  the 
following  (see  Reddy  [7]): 


s;n=o 

in  Q 

(1) 

v-Q+</  =  o 

in  Q 

(2) 

d£M  -  Q  =  0 

in  Q 

(3) 

N  -  Ad£u  -  B  dE0  =  0 

in  Q 

(4) 

M-Bd£u-DdE*D  =  0 

in  Q 

(5) 

Q- A(Vw  +  <!>)  =  0 

in  Q 

(6) 

For  free  vibration  analysis,  the  loads  are  set  to  zero  and  the  variables  assume  a  periodic  solution 
in  time,  with  a  frequency  co.  Hence,  the  governing  equations  by  the  FSDT  for  free  vibration 
analysis  of  laminated  composite  plates  are  in  turn,  as  follows  (see  Reddy  [7]): 

d£  N  +  co2I0u  +  0)2IxQ>  =  0 

in  Q 

(V) 

V  ■  Q  +  co2I0w  =  0 

in  Q 

(8) 

d£  M  -  Q  +  (o2Ij\  +  arl2<$> 

=  0  in  Q 

(9) 

N-Ad£u-Bd;;a>  =  0 

in  Q 

(10) 

M-Bd£u-Dd£<I>  =  0 

in  Q 

(ID 

Q  -  A(Vw  +  <!>)  =  0 

in  Q 

(12) 

For  convenience,  a  compact  notation  is  applied  here  which  proves  to  be  rather  useful  to  develop 
the  least-squares  functional  afterwards.  Overall,  for  static  analysis  the  governing  equations 
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include  the  plate  equilibrium  equations  in  Eqs.  (l)-(3)  and  the  laminate  linear  constitutive 
equations  in  Eqs.  (4)-(6),  while  for  free  vibration  analysis  the  governing  equations  include  the 
plate  equations  of  periodic  motion  in  Eqs.  (7)-(9)  and  the  same  laminate  constitutive  equations 
repeated  in  Eqs.  ( 1 0)-(  12),  all  written  in  tenns  of  the  independent  variables.  Specifically,  in  the 
current  notation,  the  in-plane  displacements  u,  the  transverse  deflection  w,  the  rotations  O,  the 
in-plane  force  resultants  N,  the  moment  resultants  M  and  the  transverse  force  resultants  Q,  are 
assumed  to  be  in  the  fonn  specified  below: 


u  =  [«0 

v0]r,  N  =  N„  NJt 

(13a) 

w=w0, 

m  =  [m„ 

(13b) 

®=k 

d\Q  =  \q,  0,]t 

(13c) 

In  addition,  the  fonn  of  the  differential  operator  used  repeatedly  in  the  previous  governing 
equations  is  given  by: 
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d/dx  0  d/dy 
0  d/dy  d/dx 


(14) 


In  the  laminate  constitutive  equations,  equally  written  in  Eqs.  (4)-(6)  and  in  Eqs.  ( 1 0)-(  12),  a 
matrix  form  for  the  laminate  stiffnesses  is  considered  as  follows: 


B  = 
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-m2  ^22 

_^16  ^26 
B\\  -®12 

-®12  B22 
-®16  B26 
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A  = 

K  Ass 

_^45 

^45 

A44 

Dn 

A 

d16 

D  = 

D12 

D22 

D26 

A 6 

D26 

D66 

(15a) 


(15b) 


Here,  Ay  are  called  the  extensional  stiffnesses,  Dy  the  bending  stiffnesses  and  By  the  bending- 
extensional  coupling  stiffnesses.  The  factor  K  represents  the  shear  correction  coefficient,  which 
takes  the  standard  value  5/6  in  the  later  numerical  examples.  Specifically,  the  laminate 
stiffnesses  Ay,  By  and  D y  are  in  turn  defined  in  terms  of  the  lamina  stiffnesses,  i.e.  the  lamina 
plane-stress  reduced  stiffnesses  transfonned  to  the  xv-plane  of  the  laminate,  as  shown  (see  Reddy 
[7]): 


Aij  =YJQyk\zk+\  ~zk) 

k=l 

(16a) 

1  N  .  . 

bu  =^Y.Qlk)(zM-zi) 

Z  k= 1 

(16b) 

1  N  ,  . 

D  k= 1 

(16c) 
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Furthermore,  the  mass  moments  of  inertia  /,  introduced  in  the  plate  equations  of  periodic  motion 
in  Eqs.  (7)-(9)  are  defined  in  terms  of  the  laminas  density,  as  also  shown: 


h  =Yjp'i\zk+x-zk) 

k=\ 

(17a) 

a  =\npt](zl+i  -zl) 

^  k=i 

(17b) 

i  n  ,  , 

/2=(Zrf)(4,-^5) 

^  k= 1 

(17c) 

For  completeness,  the  proper  boundary  conditions  for  all  possible  support  types  used  in  a 
rectangular  laminated  composite  plate,  as  illustrated  in  Fig.  1,  are  now  specified: 


x  =  0,a : 

=  %  =  <Py  =  =  Mxx  =  0 

on  rsl 

(18a) 

U0  =w0=  <py  =  Nxy  =  Mxx  =  0 
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(18b) 
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on  rF 

(18c) 
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Clearly,  combinations  of  the  boundary  conditions  in  Eqs.  (18)  and  (19)  can  be  made,  in  view  of 
the  support  types  considered  for  each  of  the  four  edges  of  the  rectangular  laminate.  The  notation 
used  hereafter  for  the  boundary  conditions  is  such  that  each  edge  is  specified  as  simply  supported 
(S),  free  (F)  or  clamped  (C),  strictly  in  this  sequence:  x  =  0 ,  x  =  a  ,  y  =  0  and  y  =  b.  For  simply 
supported  boundary  conditions  two  types  of  support  are  possible,  usually  named  SI  and  S2. 
Whichever  is  being  used  is  specified  in  the  end.  Hence,  the  notation  FFSS 1  for  example,  is  used 
to  denote  a  rectangular  laminate  for  which  the  edges  x  =  0,a  are  free  and  the  edges  y  =  0,b  are 
simply  supported  of  type  1 . 

2.3  Least-Squares  Formulation 

From  a  practical  standpoint,  it  is  best  to  develop  a  least-squares  finite  element  model  that  allows 
the  use  of  C°  basis  functions  in  order  to  reduce  the  higher  regularity  requirements  of  any 
weighted  residual  fonnulation  (see  Pontaza  [9]).  Therefore,  whenever  necessary  the  governing 
equations  should  be  transformed  into  an  equivalent  first-order  system,  which  implies  that 
additional  independent  variables  need  to  be  introduced.  Nonetheless,  this  transformation  can  be 
argued  to  be  somewhat  beneficial,  as  the  auxiliary  variables  may  represent  physically  meaningful 
variables,  in  the  framework  of  mixed  formulations.  In  the  present  case  though,  both  systems  of 
governing  equations  are  already  of  first-order,  namely  the  system  given  by  Eqs.  ( 1  )-(6)  for  static 
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analysis  and  the  system  given  by  Eqs.  (7)-(  12)  for  free  vibration  analysis.  Hence,  it  is  only 
necessary  to  develop  the  least-squares  functional  appropriate  to  each  analysis  and  minimize  it 
with  respect  to  the  chosen  approximation  spaces  to  obtain  the  correspondent  least-squares  finite 
element  model. 

Basically,  the  least-squares  functional  is  defined  by  measuring  the  residuals  of  the 
governing  equations  in  terms  of  suitable  norms.  To  do  so,  standard  notation  is  used.  Specifically, 
the  norm  corresponding  to  the  Sobolev  space  s>0  is  denoted  by  ||sD  and  H'(o) 

represents  the  product  space  [hs  (fl)]" ,  where  n  is  the  number  of  space  dimensions. 


Thus,  the  least-squares  functional  for  the  static  analysis  of  laminated  composite  plates  is 
based  on  the  norms  of  Eqs.  (l)-(6),  as  follows: 
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(20) 


Similarly,  the  least-squares  functional  for  free  vibration  analysis  is,  in  turn,  based  on  the  norms 
of  Eqs.  (7)-(12),  as  follows: 


Jv  (u,w,0>,N,M,Q;ry)  =  \  +  ct)2I0  u  +  ry2/,® 
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Accordingly,  the  least-squares  principle  can  be  stated  as: 

Find  (u,w,<D,N,M,Q)e  X  such  as  for  all  (s,t,'F,0,P,R)e  X 


(21) 


Js(u,  w,®,N,M,Q;g)  <  if  static  analysis  (22) 

J v  (u,  w,<D,N,M,Q;tfj)  <  «/K(s,t,'P,0,P,R;  (O ) ,  if  free  vibration  analysis  (23) 


The  space  X is  defined  below  and  satisfies  the  support  type  boundary  conditions: 

X  =  {(u,  w,®,N,M,Q)g  H'(fi)xE1(fi)xH'(Q)xHl(fi)xH1(o)xH1(Q)}  (24) 

Hence,  the  least-squares  formulation  leads  the  static  and  free  vibration  analysis  of  laminated 
composite  plates  to  the  unconstrained  minimization  problems  given  by  Eq.  (22)  and  Eq.  (23), 
respectively.  Subsequently,  the  finite  element  models  are  developed  by  minimizing  the  least- 
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squares  functional.  Specifically,  the  Euler-Lagrange  equations  are  derived  for  each  minimization 
problem,  so  as  to  obtain  the  least-squares  variational  problem  for  static  and  free  vibration 
analysis. 

2.4  Finite  Element  Models 

The  mentioned  least-squares  variational  problems  give  rise  to  the  corresponding  finite  element 
models  for  static  and  free  vibration  analysis  of  laminated  composite  plates.  Accordingly,  the 
infinite  dimensional  space  X  is  now  restricted  to  the  finite-dimensional  subspace  X/,p,  where  h 
denotes  the  mesh  parameter  and  p  the  order  for  the  variables  basis  functions. 

Ultimately,  the  mixed  least-squares  finite  element  model  for  static  analysis  takes  the  following 
matrix  form: 
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In  addition,  the  mixed  least-squares  model  for  free  vibration  analysis  develops  into  a  quadratic 
eigenvalue  problem  as  follows,  also  in  matrix  form: 
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Both  finite  element  models  yield  only  symmetric  matrices  by  means  of  the  least-squares 
formulation.  Specifically,  in  view  of  the  adopted  FSDT  mixed  formulation,  all  matrices  can  be 
structured  in  13x13  submatrices  by  considering  the  variables  separately:  namely,  the  5 
generalized  displacements  and  the  8  stress  resultants.  The  explicit  integral  expressions  of  all 
nonzero  submatrices  Kip  C(/-  and  My  and  all  nonzero  subvectors  F,  are  included  in  Appendix  A. 
The  stiffness  matrix  K  is  shared  by  both  models  and  it  is  not  only  symmetric  but  also  positive- 
definite,  once  the  boundary  conditions  are  properly  imposed.  This  fact  allows  the  use  of  robust 
solvers  for  the  static  analysis  of  laminated  composite  plates  by  the  mixed  least-squares  model.  In 
the  quadratic  eigenvalue  problem,  besides  the  stiffness  matrix  K,  both  C  and  M  appear  as  indeed 
mass  matrices.  The  difference  is  that  the  matrix  M  refers  to  the  mass  relation  among  generalized 
displacements  only,  whereas  the  matrix  C  translates  the  mass  coupling  between  generalized 
displacements  and  stress  resultants. 

In  view  of  the  finite  element  method,  the  approach  for  numerically  evaluating  the 
integrals  implicit  in  Eqs.  (25)  and  (26)  is  to  map  the  finite  elements  that  form  the  entire  model 
into  a  master  element.  For  both  least-squares  models  implemented,  the  integrals  are  evaluated 
using  full  integration  through  Gauss  quadrature  rules,  which  implies  a  master  element  with  the 
coordinate  system:  —  1  <  (^,  77)  <  1 .  Over  this  master  element,  the  variables  are  approximated  by 
equal  high-order  C°  basis  functions  as  exemplified  for  the  transverse  deflection  below: 
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(27) 


%  i^n)  «  w*p  (£  77)  =  X  wj<pJ  (£  *7) 

M 

Here,  wj  denotes  nodal  values  for  the  transverse  deflection,  (p/  the  associated  high-order  C°  basis 
functions  and  n  the  number  of  nodes  in  the  finite  element. 

Specifically,  the  mixed  least-squares  model  for  static  analysis  uses  the  classical  C° 
Lagrange  basis  functions,  whereas  the  later  model  for  free  vibration  analysis  uses  instead  C° 
interpolant  polynomials  of  Gauss-Lobatto-Legendre  quadrature  points.  The  last  basis  functions 
were  initially  used  in  the  spectral  element  method  and  are  in  fact  more  suitable  for  high-order 
finite  elements.  In  any  case,  the  two-dimensional  basis  functions  are  given  by  tensor  products  of 
the  corresponding  one-dimensional  basis  functions. 

The  well-known  one-dimensional  Lagrange  basis  functions  of  order  p  =  N- 1  can  be  defined  by  N 
equally  spaced  nodes  ^  ,  given  =  -1  and  %N  =  1 ,  as  follows: 

<28> 

7=1  G  bj 
j*i 


Alternatively,  the  one-dimensional  basis  functions  derived  from  Gauss-Lobatto-Legendre  points 
of  order  p  =  N- 1,  can  be  written  using  the  Legendre  polynomial  of  same  order  PN.\,  as  follows: 


w(f)= 


Ayv-iK-,  Of-#,) 


(29) 


where,  represents  N nodes  now  given  by  the  roots  of  (l-^2)/^_,(^)=0  in  the  interval  [-1,1]. 
For  more  details  on  these  basis  functions  see  Warburton  et  al  [14]. 

Furthermore,  the  implemented  least-squares  models  for  static  and  free  vibration  analysis  differ 
also  in  the  Gauss  quadrature  rule  used.  Gauss-Legendre  rule  is  used  for  static  analysis,  whereas 
Gauss-Lobatto-Legendre  rule  is  conveniently  employed  for  free  vibration  analysis,  due  to  the 
chosen  basis  functions. 

The  global  system  of  equations  either  for  static  or  free  vibration  analysis  is  then 
assembled  from  the  element  contributions  by  the  standard  summation  approach,  followed  by 
imposition  of  the  appropriate  boundary  conditions  (see  Reddy  [6]).  In  fact,  unlike  weak  form 
finite  element  models  that  allow  weak  imposition  of  stress  resultants  by  integral  boundary  terms, 
the  least-squares  models  only  allow  strong  imposition  of  the  boundary  conditions  both  for  stress 
resultants  and  generalized  displacements. 

2.5  Computational  Specifics 

The  mixed  least-squares  model  for  static  analysis  of  laminated  composite  plates  yields  a 
symmetric  and  positive-definite  system  of  algebraic  equations.  Hence,  its  numerical  solution  is  in 
fact  straightforward.  However,  post-computation  of  strains  and  stresses  such  that  no  numerical 
differentiation  is  carried  out  requires  a  more  attentive  procedure. 

In  detail,  the  membrane  and  flexural  strains  are  computed  first  using  the  following 
laminate  constitutive  relations,  once  the  solution  of  stress  resultants  is  known: 
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B  D 


Q  =  Ae° 


In  these  equations,  the  implied  vector  form  for  the  stress  resultants  is  previously  defined  in  Eq. 
(13)  and  the  matrix  fonn  for  the  laminate  stiffnesses  is  given  by  Eq.  (15).  Hence,  only  the 
appropriate  form  for  the  strain  components  needs  to  be  specified,  as  follows: 


Secondly,  the  in-plane  and  transverse  stresses  are  computed  given  the  prior  membrane  and 
flexural  strains,  again  through  the  laminate  constitutive  relations  but  yet  in  another  form: 

=qW(£°  +z£1)?  =Q(iV  (32) 


Specifically,  the  in-plane  stresses  are  computed  on  the  top  and  bottom  of  each  Mi  layer  while  the 
transverse  stresses  are  only  computed  within  each  Mi  layer,  in  agreement  with  the  FSDT  stress 
variations  through  the  laminate  thickness.  For  clearness,  in  the  previous  equations,  it  is 
considered  the  following  matrix  form  for  the  lamina  stiffnesses: 
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Additionally,  the  stress  components  are  defined  in  a  similar  manner  as  the  strains  components 
before,  as  follows: 


CTW=[<7W  <7W  <7wr,  6W=[<7(t)  <7W1 

L  xx  yy  xy  \  ’  L  xz  yz  J 


(34) 


The  post-computation  described  for  static  analysis  of  laminated  composite  plates  ensures  that  the 
computed  stresses  experience  no  loss  of  accuracy  through  differentiation  and  are  evaluated  in 
nodal  points  as  any  other  variable. 

Regarding  the  free  vibration  analysis  of  laminated  composite  plates,  the  mixed  least- 
squares  model  yields  a  quadratic  eigenvalue  problem  involving  symmetric  matrices.  Since 
numerical  algorithm  design  for  quadratic  eigenproblems  is  still  an  active  research  topic,  the  main 
endeavor  is  to  pursue  an  efficient  method  to  solve  the  quadratic  eigenproblem  under 
consideration. 

One  approach  is  to  use  methods  that  tackle  the  quadratic  eigenproblem  directly,  usually 
variants  of  Newton’s  method  that  find  one  eigenpair  at  a  time.  However,  the  availability  of  such 
methods  is  rather  deficient.  The  approach  actually  chosen  is  to  transform  the  quadratic 
eigenproblem  into  an  equivalent  “linear”  generalized  eigenproblem  to  allow  the  use  of  traditional 
methods  for  solution  of  eigenvalue  problems  (see  Bai  et  al  [15]).  These  methods  are  in  fact 
widely  available.  In  detail,  the  current  implementation  uses  ARPACK  subroutines,  which  have  a 
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long  proven  robustness  and  accuracy,  to  compute  a  few  eigenvalues  and  eigenvectors  with 
Implicitly  Restarted  Amoldi  Methods  (IRAM). 

Specifically,  the  assembled  quadratic  eigenvalue  problem  is  of  the  following  form: 


(M+T[c]+T2[m]){a}={0},  A  =  co2  (35) 

For  reasons  soon  made  clear,  an  invert  spectral  transformation  is  first  considered,  as  follows: 
([M]  +  //[c]  +  //2[r]){A}={0},  jU  =  \/A  (36) 


Then,  the  desired  transformation  into  an  equivalent  generalized  eigenproblem  takes  the  form 
specified  below: 

[A]{x}=  h[b]{x]  (37) 

where, 


U] 


0 

-M 


I 

0 


(38) 


Ultimately,  this  approach  reduces  the  original  quadratic  eigenproblem  into  a  non-symmetric 
generalized  eigenproblem,  where  the  matrix  B  is  still  symmetric  and  positive-definite.  In  fact,  the 
reason  for  the  prior  invert  transformation  is  to  ensure  that  the  matrix  B  is  positive-definite  by 
depending  on  the  stiffness  matrix  K  rather  than  the  mass  matrix  M  (besides  the  identity  matrix  I). 
Considering  the  above  transformations,  the  equivalent  generalized  eigenproblem  can  be 
efficiently  solved  by  ARPACK  subroutines.  In  addition,  the  few  computed  eigenvalues  and 
eigenvectors  of  the  equivalent  eigenproblem  must  then  be  used  to  obtain  the  eigenvalues  and 
eigenvectors  of  the  original  quadratic  eigenproblem  through  the  relations  above. 

Finally,  the  specific  nature  of  the  finite  element  matrices  K,  C  and  M,  render  the  quadratic 
eigenproblem  and  the  equivalent  generalized  eigenproblem  generally  complex  solutions. 
However,  once  the  original  eigenvalues  and  eigenvectors  are  recovered,  the  complex  solutions 
(in  conjugate  pairs)  show  a  negligible  imaginary  part  relatively  to  the  real  counterpart.  So,  for 
practical  purposes,  only  the  real  part  of  the  solution  of  free  vibration  analysis  of  laminated 
composite  plates  is  reported  in  the  following  numerical  examples. 

2.6  Numerical  Examples 

The  predictive  capabilities  of  the  proposed  mixed  least-squares  models  are  now  demonstrated 
through  selected  problems  of  static  and  free  vibration  analysis  of  laminated  composite  plates. 
Specifically,  four  square  laminated  composite  plates  are  considered  with  different  boundary 
conditions  and  a  range  of  side-to-thickness  ratios,  covering  thick  to  thin  laminates.  The  particular 
laminates  under  analysis  are  two  cross-ply  laminates  (0/90)  and  (0/90/0/90/0)  and  two  angle-ply 
laminates  (-45/45)  and  (30/-60/60/-30).  In  addition,  the  selected  problems  of  static  analysis 
concern  laminated  composite  plates  under  a  unifonnly  distributed  load  of  intensity  q0 . 

The  material  properties  for  all  layers  of  the  given  laminates  are  assumed  to  be  the  same, 
as  shown: 
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El  /E2  =  25 ,  Gn  =  Gu  =  0.5 E2 ,  G2,  =  0.2 E2 ,  v12  =  0.25  (39) 

Moreover,  the  subsequent  numerical  results  both  in  graphical  and  tabular  forms  for  the  main 
variables  are  nondimensionalized,  as  specified  below: 
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Analytical  solutions  are  also  presented  alongside  the  numerical  results  for  comparison.  Basically, 
the  well-known  FSDT  Navier  solutions  or  Levy  solutions  for  static  or  free  vibration  analysis  of 
laminated  composite  plates  are  reported,  according  to  the  problem  under  analysis  (see  Reddy 
[7]).  For  static  analysis,  the  Navier  series  are  evaluated  for  m,«  =  l,...,40  and  the  Levy  series  for 
«  =  1,...,40 . 


Cross-ply  laminates 

The  first  selected  problem  is  the  static  analysis  of  the  antisymmetric  laminate  (0/90)  with  SSSS1 
boundary  conditions.  A  unifonn  mesh  of  4x4  square  elements  is  used  to  model  the  composite 
plate  and  an  increasing  order  for  the  mixed  least-squares  element  is  considered.  Specifically,  4th, 
6th  and  8th-order  elements  are  successively  applied  in  order  to  investigate  the  ^-convergence  of 
the  proposed  model.  The  computed  results  are  summarized  in  Table  1. 


Table  1.  Static  results  for  the  laminate  (0/90)  SSSS1  using  a  uniform  mesh  4x4. 


a/h 

/7-order 

M  (-  -) 

Mrv(0,0) 

Qyib  0) 

^ yy  \2  5  2  / 

<>(0,0) 

10 

4 

1.9469 

0.6265 

-0.1621 

3.4675 

1.0712 

0.0973 

0.5944 

6 

1.9469 

0.6268 

-0.1605 

3.4698 

1.0716 

0.0963 

0.5948 

8 

1.9469 

0.6268 

-0.1604 

3.4703 

1.0716 

0.0962 

0.5949 

Analytical 

1.9469 

0.6268 

-0.1603 

3.4194 

1.0716 

0.0962 

0.5862 

20 

4 

1.7583 

0.6288 

-0.1609 

3.4801 

1.0745 

0.0965 

0.5966 

6 

1.7582 

0.6291 

-0.1580 

3.4883 

1.0748 

0.0948 

0.5980 

8 

1.7582 

0.6291 

-0.1576 

3.4880 

1.0748 

0.0946 

0.5979 

Analytical 

1.7582 

0.6290 

-0.1575 

3.4369 

1.0748 

0.0945 

0.5892 

100 

4 

1.6980 

0.6297 

-0.1621 

3.4560 

1.0757 

0.0973 

0.5924 

6 

1.6980 

0.6301 

-0.1571 

3.4938 

1.0762 

0.0943 

0.5989 

8 

1.6981 

0.6301 

-0.1561 

3.4936 

1.0763 

0.0937 

0.5989 

Analytical 

1.6980 

0.6300 

-0.1558 

3.4427 

1.0762 

0.0935 

0.5902 
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The  effect  of  mesh  refinement  is  examined  as  well  using  only  4th-order  elements  and  a  uniform 
mesh  of  4x4,  5x5  and  8x8  square  elements.  These  numerical  results  are  in  turn  shown  in  Table  2. 
Both  tables  include  the  results  for  transverse  deflection,  stress  resultants  and  stresses,  considering 
side-to-thickness  ratios  of  10,  20  and  100.  In  particular,  for  all  stress  results  tabulated  henceforth 
the  right  superscript  specifies  the  layer  where  the  results  are  referred  to  and  when  necessary 
also  the  bottom  or  top  interface  by  b  or  t,  respectively.  Furthermore,  Navier  analytical  solutions 
for  this  problem  are  shown  throughout  both  tables. 

Table  2.  Static  results  for  the  laminate  (0/90)  SSSS1  using  the  4th-order  element. 


a/h 

Mesh 

M  (-  -) 

yy\2>2/ 

Mj 0,0) 

Qy{ f.0) 

a(2nU 

yy  \2,2> 

^(0,0) 

10 

4x4 

1.9469 

0.6265 

-0.1621 

3.4675 

1.0712 

0.0973 

0.5944 

5x5 

1.9469 

0.6268 

-0.1612 

3.4709 

1.0716 

0.0967 

0.5950 

8x8 

1.9469 

0.6267 

-0.1605 

3.4707 

1.0716 

0.0963 

0.5950 

Analytical 

1.9469 

0.6268 

-0.1603 

3.4194 

1.0716 

0.0962 

0.5862 

20 

4x4 

1.7583 

0.6288 

-0.1609 

3.4801 

1.0745 

0.0965 

0.5966 

5x5 

1.7583 

0.6291 

-0.1594 

3.4886 

1.0748 

0.0957 

0.5980 

8x8 

1.7582 

0.6290 

-0.1580 

3.4882 

1.0748 

0.0948 

0.5980 

Analytical 

1.7582 

0.6290 

-0.1575 

3.4369 

1.0748 

0.0945 

0.5892 

100 

4x4 

1.6980 

0.6297 

-0.1621 

3.4560 

1.0757 

0.0973 

0.5924 

5x5 

1.6980 

0.6300 

-0.1598 

3.4968 

1.0762 

0.0959 

0.5995 

8x8 

1.6980 

0.6300 

-0.1572 

3.4926 

1.0762 

0.0943 

0.5987 

Analytical 

1.6980 

0.6300 

-0.1558 

3.4427 

1.0762 

0.0935 

0.5902 

Overall,  the  numerical  results  in  Tables  1  and  2  are  in  good  agreement  with  the  analytical 
solutions  for  the  entire  range  of  side-to-thickness  ratios  analyzed.  In  fact,  the  centre  transverse 
deflection  is  rightly  predicted  using  just  4th-order  elements  even  when  a  thin  laminate  is 
considered  (with  a  side-to-thickness  ratio  of  100).  Most  notably,  convergence  towards  the 
analytical  solution  is  verified,  either  with  p-  or  //-refinements,  as  should  be  expected.  Actually, 
an  exact  study  on  the  asymptotic  behaviour  of  mixed  finite  elements  based  on  least-squares 
formulation  can  be  found  in  Pontaza  [9],  In  the  case  of  the  transverse  force  resultant  and 
transverse  stress,  it  is  noted  a  slight  discrepancy  between  the  analytical  and  numerical  results, 
which  is  believed  to  be  due  to  a  not  as  accurate  analytical  representation. 

In  Fig.  2,  it  is  shown  the  predictive  distributions  of  the  in-plane  stresses  ovy  and  axy  along 
the  laminate  thickness,  using  a  side-to-thickness  ratio  of  10.  The  proper  Navier  solution  is  plotted 
alongside  the  numerical  results  obtained  with  the  8th-order  element  in  a  4x4  mesh. 
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Figure  2.  In-plane  stresses  for  the  laminate  (0/90)  SSSS1  with  a/h  =  10  . 

Continuing  the  analysis  of  the  same  antisymmetric  laminate  (0/90)  with  SSSS1  boundary 
conditions,  the  problem  of  free  vibration  is  considered  next.  In  this  case,  a  uniform  mesh  of  4x4 
square  elements  is  used  together  with  the  4th-order  mixed  least-squares  element.  The  results  for 
the  fundamental  frequency  using  the  same  side-to-thickness  ratios  as  before  are  presented  in 
Table  3,  along  with  the  Navier  solutions. 

Table  3.  Free  vibration  results  for  the  laminate  (0/90)  SSSS1 
using  a  uniform  mesh  4x4. 


a/h 

/7-order 

cox 

10 

4 

8.9006 

Analytical 

8.9001 

20 

4 

9.4746 

Analytical 

9.4745 

100 

4 

9.6873 

Analytical 

9.6873 

As  mentioned  earlier,  the  mixed  least-squares  model  for  free  vibration  analysis  yields  a 
quadratic  eigenproblem,  that  to  be  later  solved  by  ARPACK  subroutines  requires  an  invert 
spectral  transformation  followed  by  another  transfonnation  into  an  equivalent  generalized 
eigenproblem.  With  this  mind,  the  smaller  eigenvalues  of  the  original  quadratic  eigenproblem 
correspond  to  the  larger  eigenvalues  of  the  generalized  eigenproblem.  Hence,  to  compute  the 
lower  natural  frequencies  of  the  quadratic  eigenproblem,  it  is  specified  to  ARPACK  to  extract 
the  eigenvalues  of  largest  real  part.  Specifically,  the  results  in  Table  3  are  obtained  by  specifying 
the  extraction  of  2  eigenvalues  (which  come  as  a  complex  conjugate  pair)  using  50  Arnoldi  basis 
vectors.  In  fact,  the  choice  of  both  these  numbers  for  an  optimal  performance  of  ARPACK  is 
truly  problem  dependent.  Since  this  optimum  choice  was  not  investigated,  whenever  free 
vibration  results  are  presented  the  number  of  eigenvalues  requested  and  number  of  Arnoldi  basis 
vectors  used  are  explicitly  stated. 

Returning  to  Table  3,  the  computed  fundamental  frequencies  are  remarkably  in 
agreement  with  the  analytical  solutions,  for  all  side-to-thickness  ratios  considered.  In  fact,  the 
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4th-order  mixed  least-squares  element  appears  to  be  quite  sufficient  to  obtain  very  goods  results 
for  the  fundamental  frequency  in  this  case. 

The  next  problem  is  the  static  analysis  of  the  symmetric  laminate  (0/90/0/90/0)  with 
FFSS1  boundary  conditions.  Once  more,  /^-convergence  of  the  proposed  mixed  least-squares 
element  is  inspected,  using  in  turn  4th,  6th  and  8th-order  elements  in  a  uniform  mesh  of  4x4 
square  elements.  Table  4  shows  the  results  for  the  transverse  deflection,  stress  resultants  and 
stresses,  using  again  side-to-thickness  ratios  of  10,  20  and  100.  For  this  problem,  Levy  analytical 
solutions  are  presented  for  comparison. 


Table  4.  Static  results  for  the  laminate  (0/90/0/90/0)  FFSS1  using  a  uniform  mesh  4x4. 


ajh 

/7-order 

1V1  yy  \2  ’  2  / 

aif.o) 

CF<4,)( 4.  2-) 

^  XX  \2  9  2  / 

V2  ’  2 / 

^4)(f,o) 

10 

4 

3.0600 

0.0061 

1.2458 

4.9787 

0.0179 

1.8719 

0.9335 

6 

3.0600 

0.0061 

1.2458 

4.9785 

0.0179 

1.8718 

0.9335 

8 

3.0600 

0.0061 

1.2458 

4.9785 

0.0179 

1.8718 

0.9335 

Analytical 

3.0600 

0.0061 

1.2458 

4.9279 

0.0179 

1.8718 

0.9240 

20 

4 

2.7082 

0.0069 

1.2449 

4.9770 

0.0179 

1.8705 

0.9332 

6 

2.7082 

0.0070 

1.2449 

4.9750 

0.0179 

1.8706 

0.9328 

8 

2.7082 

0.0069 

1.2449 

4.9748 

0.0179 

1.8705 

0.9336 

Analytical 

2.7082 

0.0069 

1.2449 

4.9241 

0.0179 

1.8705 

0.9233 

100 

4 

2.5956 

0.0072 

1.2446 

4.9787 

0.0179 

1.8700 

0.9335 

6 

2.5957 

0.0073 

1.2446 

4.9792 

0.0179 

1.8701 

0.9336 

8 

2.5955 

0.0074 

1.2446 

4.9783 

0.0179 

1.8701 

0.9334 

Analytical 

2.5957 

0.0074 

1.2446 

4.9230 

0.0179 

1.8701 

0.9231 

Noticeably,  the  computed  numerical  results  in  Table  4  show  excellent  agreement  with  the 
analytical  solutions,  even  more  than  in  the  previous  static  problem.  Not  only  the  centre 
transverse  deflection,  but  also  stress  resultants  and  stresses  can  be  predicted  exactly  using  just 
4th-order  elements,  for  all  side-to-thickness  ratios  considered.  Basically,  it  seems  that  most  of  the 
computed  results  are  converged  using  the  4th-order  elements,  including  the  results  for  a  thin 
laminate.  Again,  only  the  transverse  force  resultant  and  transverse  stress  show  a  small 
discrepancy  between  the  analytical  and  numerical  results.  This  occurrence  is  more  or  less 
apparent  throughout  every  static  analysis  problems  here  presented.  So,  to  avoid  repetition,  it  is 
understood  henceforth  that  the  static  analytical  solutions  for  these  variables  may  not  be  fully 
converged  in  the  series  representation. 

Fig.  3  shows  the  distributions  for  the  in-plane  stresses  axx  and  ayy  through  the  laminate 
thickness,  using  a  side-to-thickness  ratio  of  10.  These  plots  contain  the  numerical  results  by  the 
8th-order  element  in  a  4x4  mesh  with  the  Levy  solution. 
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Figure  3.  In-plane  stresses  for  the  laminate  (0/90/0/90/0)  FFSS1  with  a/h  =  10  . 

To  conclude  the  static  analysis  of  both  cross-ply  laminates  previously  considered,  the 
transverse  deflections  along  the  line  x  =  al  2  of  the  plate  are  plotted  in  Fig.  4.  Namely,  the 
antisymmetric  laminate  (0/90)  SSSS1  and  the  symmetric  laminate  (0/90/0/90/0)  FFSS1.  Fig.  4 
contains  then,  the  corresponding  analytical  solutions  and  numerical  results  using  the  6th-order 
element  in  a  4x4  mesh,  given  side-to-thickness  ratios  of  10  and  100  (i.e.  for  thick  and  thin 
laminates,  respectively). 


0.0  0.2  0.4  0.6  0.8  1.0 


y/a 

Figure  4.  Transverse  deflections  along  x  =  a/2  for  the  cross-ply  laminates. 

Evidently,  Fig.  4  demonstrates  that  the  computed  results  by  the  proposed  mixed  least- 
squares  element  for  static  analysis  are  in  any  case  well  in  agreement  with  the  analytical  solutions, 
whether  thin  or  thick  laminates  are  concerned.  In  fact,  results  so  far  suggest  that  this  mixed  least- 
squares  model  is  insensitive  to  shear  locking,  at  least  for  the  /> levels  considered. 

The  last  selected  problem  of  this  section  devoted  to  cross-ply  laminates  is  the  free  vibration 
analysis  of  the  former  symmetric  laminate  (0/90/0/90/0),  but  with  SSSS2  boundary  conditions 
instead.  This  time,  the  5  lowest  natural  frequencies  are  investigated  using  a  unifonn  mesh  of  4x4 
square  elements  and  the  4th-order  mixed  least-squares  element.  The  numerical  results  for  these 
natural  frequencies  are  given  in  Table  5  for  the  usual  side-to-thickness  ratios,  along  with  the 
Navier  solutions.  In  particular,  this  computation  uses  ARPACK  to  extract  of  a  cluster  of  14 
eigenvalues  (in  complex  conjugate  pairs)  with  57  Arnoldi  basis  vectors. 
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Table  5.  Free  vibration  results  for  the  laminate  (0/90/0/90/0)  SSSS2  using  a  unifonn  mesh  4x4. 


a/h 

/;-order 

cox 

( o2 

CO  3 

CO  ^  CO  ^ 

10 

4 

12.5650 

24.0078 

30.0457 

36.5359 

40.1397 

Analytical 

12.5651 

24.0120 

30.0549 

36.5559 

40.1329 

20 

4 

14.3875 

29.2513 

42.4381 

50.1919 

54.8037 

Analytical 

14.3875 

29.2508 

42.4421 

50.2606 

54.8754 

100 

4 

15.1909 

31.8771 

51.7589 

60.3319 

64.9051 

Analytical 

15.1909 

31.8770 

51.7581 

60.3285 

64.9772 

Overall,  the  5  computed  natural  frequencies  exhibit  an  outstanding  accordance  with  the 
analytical  solutions,  regardless  of  the  side-to-thickness  ratios  considered.  Actually,  the  4th-order 
element  is  able  to  predict  quite  well  the  natural  frequencies  higher  than  the  fundamental 
frequency,  even  though  the  results  slightly  worsen  as  the  natural  frequencies  increase  (which  is 
expected  to  some  extent). 

The  following  surface  plots  illustrate  the  modes  of  vibration  computed  for  these  5  natural 
frequencies,  concerning  only  the  transverse  deflection. 


Figure  5.  Modes  of  vibration  1  and  2  for  the  laminate  (0/90/0/90/0)  SSSS2. 
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Figure  6.  Modes  of  vibration  3  and  4  for  the  laminate  (0/90/0/90/0)  SSSS2. 
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Figure  7.  Mode  of  vibration  5  for  the  laminate  (0/90/0/90/0)  SSSS2. 

The  represented  inodes  of  vibration  in  Figs.  5-7  are  also  obtained  using  the  4th-order  mixed 
least-squares  element  in  a  4x4  uniform  mesh  and  specifically  for  a  side-to-thickness  ratio  of  10, 
although  qualitatively  the  modes  are  the  same  for  all  side-to-thickness  ratios.  To  be  precise,  all 
these  graphics  are  constructed  from  the  eigenvectors  solution  (for  transverse  deflection)  and 
taking  into  account  the  4th-order  basis  functions  as  well  as  the  finite  elements  geometry 
transformations  between  the  master  element.  So,  ultimately,  the  overall  mode  shape  is  displayed 
by  putting  together  every  finite  element  contribution. 

The  modes  of  vibration  for  the  transverse  deflection  are  indeed  in  agreement  with  the 
Navier  solutions,  as  the  specific  pair  of  harmonics  that  correspond  to  each  analytical  natural 
frequency  matches  the  represented  modes.  Explicitly,  the  pair  of  harmonics  along  x  and  y 
respectively,  for  increasing  natural  frequencies  are  (1,1),  (1,2),  (2,1),  (2,2)  and  (1,3).  Therefore, 
the  proposed  mixed  least-squares  model  for  free  vibration  analysis  based  on  C°  high-order  basis 
functions,  is  in  fact  capable  of  good  predictions  for  the  natural  frequencies  as  well  as  the  modes 
of  vibration. 

Angle-ply  laminates 

This  section  starts  with  the  static  analysis  problem  of  the  antisymmetric  laminate  (-45/45)  with 
SSSS2  boundary  conditions.  Similarly  to  the  first  selected  problem,  the  effect  of  both  p-  and  h- 
refinements  is  examined.  Table  6  shows  the  computed  results  using  a  fixed  unifonn  mesh  of  4x4 
square  elements  and  an  increasing  order  for  the  mixed  least-squares  element,  whereas  Table  7 
presents  the  results  using  a  fixed  /7-lcvcl  with  the  4th-order  element  and  increasingly  refined 
meshes.  Both  tables  include  results  of  transverse  deflection,  stress  resultants  and  stresses,  with 
side-to-thickness  ratios  of  10,  20  and  100.  The  appropriate  Navier  analytical  solutions  are 
presented  as  well. 
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Table  6.  Static  results  for  the  laminate  (-45/45)  SSSS2  using  a  uniform  mesh  4x4. 


a/h 

/7-order 

M  (-  -) 

yy\  2  ’  2  / 

Mxy{0,0)  Qy{ f,0)  afXX 

’  2/  u  xy 

(0,0)  <'( f  ,o) 

10 

4 

1.2792 

0.3718 

-0.4417 

3.2730 

0.3476 

0.4311 

0.3928 

6 

1.2792 

0.3720 

-0.4388 

3.2778 

0.3477 

0.4285 

0.3933 

8 

1.2792 

0.3720 

-0.4383 

3.2785 

0.3478 

0.4281 

0.3934 

Analytical 

1.2792 

0.3720 

-0.4379 

3.2276 

0.3477 

0.4277 

0.3873 

20 

4 

1.0907 

0.3744 

-0.4500 

3.2468 

0.3496 

0.4380 

0.3896 

6 

1.0907 

0.3745 

-0.4487 

3.2507 

0.3497 

0.4370 

0.3901 

8 

1.0907 

0.3745 

-0.4482 

3.2514 

0.3497 

0.4365 

0.3902 

Analytical 

1.0907 

0.3744 

-0.4477 

3.2005 

0.3497 

0.4360 

0.3841 

100 

4 

1.0305 

0.3762 

-0.4547 

3.2357 

0.3512 

0.4419 

0.3883 

6 

1.0305 

0.3755 

-0.4558 

3.2392 

0.3505 

0.4430 

0.3887 

8 

1.0305 

0.3743 

-0.4556 

3.2396 

0.3505 

0.4428 

0.3887 

Analytical 

1.0305 

0.3755 

-0.4549 

3.1884 

0.3505 

0.4422 

0.3826 

Table  7.  Static  results  for  the  laminate  (-45/45)  SSSS2  using  the  4th-order  element. 


a/h 

Mesh 

M  (-  -) 

yy\  2  ’  2  / 

4^  (0,0)  QX f,0)  o%% 

«)  a[lb) l 

(0,0)  <7«(f  ,o) 

10 

4x4 

1.2792 

0.3718 

-0.4417 

3.2730 

0.3476 

0.4311 

0.3928 

5x5 

1.2792 

0.3721 

-0.4405 

3.2798 

0.3478 

0.4301 

0.3936 

8x8 

1.2792 

0.3720 

-0.4391 

3.2788 

0.3477 

0.4288 

0.3935 

Analytical 

1.2792 

0.3720 

-0.4379 

3.2276 

0.3477 

0.4277 

0.3873 

20 

4x4 

1.0907 

0.3744 

-0.4500 

3.2468 

0.3496 

0.4380 

0.3896 

5x5 

1.0907 

0.3745 

-0.4496 

3.2513 

0.3497 

0.4377 

0.3902 

8x8 

1.0907 

0.3744 

-0.4487 

3.2516 

0.3497 

0.4370 

0.3902 

Analytical 

1.0907 

0.3744 

-0.4477 

3.2005 

0.3497 

0.4360 

0.3841 

100 

4x4 

1.0305 

0.3762 

-0.4547 

3.2357 

0.3512 

0.4419 

0.3883 

5x5 

1.0305 

0.3754 

-0.4548 

3.2355 

0.3504 

0.4421 

0.3883 

8x8 

1.0306 

0.3755 

-0.4551 

3.2393 

0.3506 

0.4424 

0.3887 

Analytical 

1.0305 

0.3755 

-0.4549 

3.1884 

0.3505 

0.4422 

0.3826 

Again,  the  numerical  results  shown  in  both  Tables  6  and  7  are  altogether  in  good 
agreement  with  the  analytical  solutions,  for  the  range  of  side-to-thickness  ratios  analyzed. 
Specifically,  convergence  of  the  numerical  results  is  once  more  verified  for  both  p-  and  h- 
refinements.  Actually,  it  is  interesting  to  note  that  the  number  of  degrees  of  freedom  given  by  the 
8th-order  elements  in  a  4x4  mesh  is  exactly  the  same  as  the  4th-order  elements  in  an  8x8  mesh. 
A  comparison  between  these  two  cases  suggests  that  /^-refinement  is  somewhat  more  efficient, 
especially  in  view  of  side-to-thickness  ratios  of  10  and  20.  However,  for  an  exact  study  on  this 
matter  see  Pontaza  [9].  Still,  in  this  problem  it  appears  that  the  4th-order  element  is  sufficient  to 
provide  accurate  predictions  for  the  transverse  deflection  as  well  as  reasonable  predictions  for 
stresses  and  stress  resultants,  whether  thin  or  thick  laminates  are  considered. 

The  in-plane  stresses  avy  and  oxy  through  the  laminate  thickness  are  plotted  in  Fig.  8, 
given  a  side-to-thickness  ratio  of  10,  and  using  the  results  obtained  by  the  8th-order  element  in  a 
4x4  mesh  together  with  the  Navier  analytical  solution. 


41 


Figure  8.  In-plane  stresses  for  the  laminate  (-45/45)  SSSS2  with  a/h  =  10  . 

The  problem  of  free  vibration  analysis  of  the  antisymmetric  laminate  (-45/45)  with 
SSSS2  boundary  is  now  considered.  The  fundamental  frequency  is  computed  by  the  4th-order 
mixed  least-squares  element  in  a  4x4  uniform  mesh,  for  the  same  side-to-thickness  ratios  as 
before.  In  this  particular  computation  with  ARPACK,  it  is  specified  the  extraction  of  2 
eigenvalues  (in  a  complex  conjugate  pair)  using  50  Arnoldi  basis  vectors.  The  computed  results 
are  presented  in  Table  8  as  well  as  the  Navier  solutions. 

Table  8.  Free  vibration  results  for  the  laminate  (-45/45)  SSSS2  using  a  uniform  mesh  4x4. 


a/h 

/;-order 

G)\ 

10 

4 

10.8951 

Analytical 

10.8951 

20 

4 

11.9327 

Analytical 

11.9329 

100 

4 

12.3408 

Analytical 

12.3408 

From  Table  8,  it  is  verified  that  the  computed  fundamental  frequencies  are  almost  entirely 
coincident  with  the  analytical  solutions,  for  these  side-to-thickness  ratios.  Hence,  concurring 
with  the  previous  free  vibration  analysis  problems,  the  mixed  least-squares  model  of  4th-order  is 
quite  capable  of  obtaining  excellent  results  for  the  fundamental  frequencies. 

The  very  last  selected  problem  considers  the  static  analysis  of  the  antisymmetric  laminate 
(30/-60/60/-30)  with  SSSS2  boundary  conditions.  In  this  case,  a  uniform  mesh  of  4x4  square 
elements  is  used  along  with  either  the  4th,  6th  or  8th-order  elements.  Table  9  shows  these 
numerical  results  for  the  transverse  deflection,  stress  resultants  and  stresses,  with  the  usual  side- 
to-thickness  ratios.  The  Navier  analytical  solutions  are  likewise  shown. 
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Table  9.  Static  results  for  the  laminate  (30/-60/60/-30)  SSSS2  using  a  uniform  mesh  4x4. 


a/h 

/;-order 

w(f,f) 

M  (-  -) 

1Y±  XX  \2  ’  2  / 

-)  <7(3,)  (- 
>  2  >yj yy  \2  ' 

10 

4 

0.9261 

0.7107 

0.4410 

3.7218 

0.4702 

0.4034 

0.5423 

6 

0.9262 

0.7125 

0.4419 

3.7172 

0.4715 

0.4041 

0.5416 

8 

0.9262 

0.7125 

0.4419 

3.7178 

0.4714 

0.4041 

0.5417 

Analytical 

0.9261 

0.7125 

0.4418 

3.6668 

0.4714 

0.4041 

0.5343 

20 

4 

0.7306 

0.7224 

0.4391 

3.7279 

0.4785 

0.3947 

0.5432 

6 

0.7307 

0.7274 

0.4413 

3.7105 

0.4821 

0.3963 

0.5407 

8 

0.7307 

0.7273 

0.4412 

3.7114 

0.4821 

0.3962 

0.5408 

Analytical 

0.7307 

0.7273 

0.4412 

3.6604 

0.4821 

0.3962 

0.5334 

100 

4 

0.6671 

0.7235 

0.4364 

3.7592 

0.4796 

0.3898 

0.5478 

6 

0.6682 

0.7326 

0.4413 

3.7074 

0.4858 

0.3937 

0.5402 

8 

0.6681 

0.7326 

0.4412 

3.7104 

0.4858 

0.3936 

0.5407 

Analytical 

0.6681 

0.7325 

0.4411 

3.6594 

0.4858 

0.3936 

0.5332 

Table  9  shows  yet  again  that  the  computed  results  are  overall  well  in  agreement  with  the 
analytical  solutions.  As  previously  stated,  it  is  evident  that  even  when  thin  laminates  are 
considered  the  proposed  mixed  least-squares  finite  element  experiences  no  shear-locking,  so  far 
as  4th-  or  higher-order  elements  are  used.  Nevertheless,  in  this  particular  problem  the  numerical 
results  for  stresses  and  stress  resultants  by  the  4th-order  element  do  not  show  as  much  accuracy 
as  in  the  previous  static  problems,  especially  for  a  side-to-thickness  ratio  of  100.  In  fact,  p- 
convergence  is  quite  apparent  when  the  subsequent  results  by  the  6th-order  elements  are 
examined. 

The  predicted  in-plane  stresses  a **  and  ayv  through  the  laminate  thickness  are  plotted  in 
the  following  Fig.  9,  for  a  side-to  thickness  ratio  of  10.  Again,  the  stress  results  are  given  by  the 
8th-order  element  in  a  4x4  unifonn  mesh  along  with  the  analytical  solution. 

Similarly  to  the  cross-ply  laminates  earlier,  the  static  analysis  of  both  angle-ply  laminates  is 
concluded  with  the  plot  of  the  transverse  deflections  along  the  line  x  =  a  /  2  of  the  plate, 
considering  side-to-thickness  ratios  of  10  and  100.  So,  it  is  shown  in  the  next  Fig.  10  the  proper 
transverse  deflections  distributions  obtained  by  the  6th-order  element  in  a  4x4  uniform  mesh 
with  the  corresponding  analytical  solutions. 
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Figure  9.  In-plane  stresses  for  the  laminate  (30/-60/60/-30)  SSSS2  with  a/h  =  10 . 
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Figure  10.  Transverse  deflections  along  x  =  a/2  for  the  angle-ply  laminates. 


It  is  apparent  in  Fig.  10  that  the  numerical  results  obtained  by  the  mixed  least-squares 
element  follow  extremely  well  the  analytical  solutions,  for  all  the  cases  presented.  Furthermore, 
no  evidence  of  shear-locking  has  been  encountered  in  any  static  and  free  vibration  analysis 
problems  considered  to  date,  intended  for  assessment  of  the  proposed  mixed  least-squares 
models. 


2.7  Concluding  Remarks 

This  report  presents  mixed  least-squares  finite  element  models  for  static  and  free  vibration 
analysis  of  laminated  composite  plates  as  a  reliable  alternative  to  the  mixed  weak  form  models. 
The  theoretical  and  computational  advantages  of  the  least-squares  variational  principle  combined 
with  mixed  formulations  are  stated  from  the  start  and  verified  for  the  proposed  models. 
Explicitly,  the  least-squares  formulation  leads  to  an  unconstrained  minimization  problem,  which 
ensures  that  no  restrictive  compatibility  conditions  are  required  among  the  mixed  finite  element 
approximation  spaces.  In  addition,  the  mixed  least-squares  discrete  models,  once  the  boundary 
conditions  are  duly  imposed,  yield  a  symmetric  and  positive-definite  stiffness  matrix. 

The  proposed  mixed  least-squares  models  for  static  and  free  vibration  analysis  of  laminated 
composite  plates  consider  the  FSDT  with  generalized  displacements  and  stress  resultants  as 
independent  variables,  using  equal-order  interpolation.  In  fact,  to  ensure  a  correct  minimization 
of  the  least-squares  functional,  high-order  C°  basis  functions  and  full  integration  are  used  to 
develop  the  discrete  finite  element  models.  Specifically,  the  model  for  static  analysis  uses  the 
classical  C°  Lagrange  basis  functions  and  the  later  model  for  free  vibration  analysis  uses  C° 
interpolant  polynomials  of  Gauss-Lobatto-Legendre  quadrature  points,  which  are  more  suitable 
basis  functions  for  high-order  elements. 

The  predictive  capabilities  of  both  mixed  least-squares  models  are  assessed  by  a  selection  of 
numerical  examples  concerning  four  laminated  composite  plates  with  different  boundary 
conditions  and  a  range  of  side-to-thickness  ratios,  from  thick  to  thin  laminates.  Overall,  the 
numerical  results  for  static  and  free  vibration  analysis  show  excellent  agreement  with  the 
analytical  solutions  for  all  problems  examined.  In  the  case  of  static  analyses,  results  for 
transverse  deflection,  stress  resultants  and  stresses  are  carefully  inspected  and  specially, 
convergence  of  the  computed  results  towards  the  analytical  solutions  is  verified  for  both  p-  and 
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/? -refinements.  In  free  vibration  analyses,  the  natural  frequencies  are  remarkably  well  predicted 
and  even  the  modes  of  vibration  are  correctly  represented.  Furthennore,  both  mixed  least-squares 
models  are  shown  to  be  insensitive  to  shear-locking  when  modeling  thin  laminates. 

2.8  Appendix  for  Part  2 
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These  integral  expressions  make  use  of  the  symmetry  in  the  laminate  stiffnesses,  namely, 
Aj  =  Ajt ,  By  =  By  and  Dy  =  Djt .  In  addition,  for  any  submatrice  relating  two  given  variables  ( a 

and  b)  the  following  relations  hold,  Kf  =  Kbb  ,  Cf  =  Chj'L  and  Mf  =  Mh‘‘ ,  which  renders 
symmetry  to  the  all  finite  element  matrices. 
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